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Abstract 

We study the single particle velocity distribution for a granular fluid of inelas- 
tic hard spheres or disks, using the Enskog-Boltzmann equation, both for the 
homogeneous cooling of a freely evolving system and for the stationary state 
of a uniformly heated system, and explicitly calculate the fourth cumulant of 
the distribution. For the undriven case, our result agrees well with computer 
simulations of Brey et al. [|J. Corrections due to non-Gaussian behavior on 
cooling rate and stationary temperature are found to be small at all inelastici- 
ties. The velocity distribution in the uniformly heated steady state exhibits a 
high energy tail ~ exp(— Ac 3 / 2 ), where c is the velocity scaled by the thermal 
velocity and A ~ l/\fe with e the inelasticity. 

I. INTRODUCTION 

Most theories for rapid granular flows are based on the assumption that the particle 
velocities are distributed according to a Gaussian or Maxwell distribution. Since granular 
particles collide inelastically, this assumption is not an obvious one, however. In fact, gran- 
ular systems are typically far from equilibrium systems in the sense that an external driving 
force is necessary to maintain a stationary or periodic state. Only in systems of nearly 
elastic particles states close to equilibrium are feasible. 

Non-Gaussian behavior in rapid granular flows has been studied in several contexts. 
Taguchi and Hayakawa observed power law behavior of the tails of the velocity dis- 
tribution in computer simulations of a bed of grains fluidized by vertical vibration, and 
were able to explain their observations. Esipov and Poschel || have solved the Boltzmann 
equation for inelastic hard spheres or disks for large velocities. For the freely evolving sys- 
tem, they found a spatially homogeneous distribution / which for large velocities decays as 
/ ~ exp(— Av/vo(t)), where v (t) is the time dependent thermal velocity and A ~ 1/e a 
constant, related to the inelasticity e = 1 — a 2 , defined in terms of the coefficient of nor- 
mal restitution a. An enhanced population for large energies was also found by Brey et 
al. 0], who, based on a BGK model kinetic equation for an undriven granular gas, found 
algebraically decaying tails with diverging velocity moments of degree > 2/e. Sela and 
Goldhirsch || have numerically obtained a perturbative solution of the Boltzmann equation 
for inelastic hard spheres to orders of 0(e), 0(ek), 0(k 2 ). The order 0(e) estimates the 
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deviation from Gaussian behavior of the homogeneous solution and contributes to the rate 
of homogeneous cooling. 

Here we will present an alternative approach for solving the Enskog-Boltzmann equation 
for homogeneous single particle distributions. We calculate explicitly the fourth moment for 
a freely evolving system of inelastic hard spheres or disks, where the distribution function 
obeys a scaling form, i.e. it is time dependent via the decaying temperature T(t) only. Our 
method follows Goldshtein and Shapiro M, who calculated the fourth moment for a freely 
evolving gas of inelastic hard spheres (d = 3), but unfortunately made an error in their 
algebra. From the fourth cumulant, we obtain the corrected cooling rate. 

Next we consider a system of inelastic hard disks or spheres, where at equal times a ran- 
dom velocity is added to the velocity of each particle, referred to as uniformly heated system 
or as random acceleration model. This idea of uniform heating allows for the existence of 
a homogeneous stationary state and was first introduced for inelastic particles on a line by 
Williams and MacKintosh [[/[]. Recently, Peng and Ohta || and Puglisi at al. || have per- 
formed computer simulations of two-dimensional and one-dimensional systems, respectively. 
For the steady state we calculate the fourth cumulant of the velocity distribution, as well 
as the corresponding stationary temperature. Moreover we will present a calculation of the 
high energy tail of the distribution function for the uniformly heated system similar to the 
one presented in Ref. || for an undriven granular gas. 

We solve the equation for the single particle velocity distribution in the homogeneous 
state by expanding it in Sonine polynomials and deriving equations for its moments. As 
the calculation of the fourth moment is already quite involved, we have not attempted to 
calculate higher moments. Comparing our calculation with the one presented in Ref. |Q, 
we note that the results for the 0(e) correction to the cooling rate are very close. Whereas 
the calculation in Ref. || provides quantitative information on the distribution function 
itself, ours provides quantitative information on its moments, and shows only qualitatively 
similar behavior for the distribution. The advantage of our method, however, is that it is 
nonperturbative in the inelasticity, i.e. the moments can be obtained for all values of the 
coefficient of restitution. 

Our starting point is the nonlinear Enskog-Boltzmann equation |TIJ for the single particle 
distribution function /(r, v, t) in a dense system of inelastic hard spheres in d dimensions. 
In the absence of external forces, the homogeneous solution f(v,t) of this equation obeys 

d t f(v 1 ,t)= X a d - 1 J dv 2 / / dff(v 12 .*){^/(vr ) *)/(vr,t)-/(vi,0/(v2,0} 

= xHf,f)- (i) 

The prime on the cr integration denotes the condition V12 • <x > 0, where er is a unit 
vector along the line of centers of the colliding spheres at contact. In direct collisions of 
inelastic hard spheres with a coefficient of normal restitution a, the initial relative velocity 
v 12 follows the inelastic reflection law \\ 2 ' & = — av i2 • The gain term in (|I|) describes the 
restituting collisions, i.e. the precollision velocities (v**,V2*) yield (vi,v 2 ) as postcollision 
ones with v* 2 ■ <x = — (l/a;)vi 2 • &. Total momentum is conserved in a binary collision, and 
consequently in direct collisions 

v* = vi - ~(1 + a)(vi 2 • tr)tr 

v 2 = v 2 + ±(l + a)(v 12 -<x)<x, (2) 
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whereas v**(a) = v*(l/a) in restituting collisions. The factor 1/a 2 in the gain term orig- 
inates from the Jacobian dv**dv 2 * = (l/a)dv 1 dv 2 and from the length of the collision 
cylinder |v* 2 • cr\dt = (l/a)|v 12 • &\dt. 

Note that for the spatially homogeneous case, the only difference between the Enskog- 
Boltzmann equation for dense systems and the Boltzmann equation for dilute systems, is 
the presence of the factor x{ n )i which is the pair correlation function at contact. It accounts 
for the increased collision frequency in dense systems, caused by excluded volume effects. 

For later reference we will also quote the equation for the rate of change of the average 
(ip) = (1/n) / dv0(v)/(v, t), where the density n = J dv/(v, t). From QTJ) it follows as 

^d? = / dVldV2 / d * (Vl2 ' &)f{vu t)/(V2 ' t] A[ ^ (Vl) + ^ (V2)] ' (3) 

where A^(vj) = ^(v*) — ^(v*) is the ip change in a direct collision. 

In the next section we study the solution of (HD for a freely evolving fluid. In the 
subsequent section a uniformly heated system of inelastic particles will be considered. 



II. HOMOGENEOUS COOLING STATE 

For the freely evolving granular fluid, Goldshtein and Shapiro || have shown that Eq. 
(H) admits an isotropic scaling solution, describing the homogeneous cooling state, with a 
single particle distribution function depending on time only through the temperature T(t) 
as 

where the thermal velocity Vo(t) is defined in terms of the temperature by T(t) = ^mv^t), 
with 

\dnT{t) = J dv\mv 2 f{y, t), (5) 

and m the particle mass. Choosing ip = ^mvf in Eq. (§) we obtain for the rate of change of 
the temperature 

^ = -^mxna^vl = -2 7 ^ T. (6) 

Here ujq is the Enskog collision frequency for elastic hard spheres, defined as the average loss 
term in Eq. ([!]), 



uj = X na d 1 



/'d<T(v 12 -<T)\ =^L X na d ^v , (7) 
J I o v 2tt 



where (. . .)o denotes an average over Maxwellian velocity distributions for vi and v 2 at 
temperature T = ^ mv o an d — 2n d / 2 /T(d/2) is the surface area of a <i-dimensional unit 
sphere. The second equality in (|6|) defines the time independent dimensionless cooling rate 
as 7 = (\f2~K / dVL d )jjL2, where 
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Hp 



dc lC ?/(/,/) 



(8) 



are the moments of the dimensionless collision integral 



1 r, 



/(/, /) = dc 2 d<x(c 12 ■ 6-) —f( c ?)f(c?) - /(d)/(ca 



(9) 



with c = v/t> (t). Using Eqs. (H) and (Q), the scaling form /(c) satisfies the integral equation 

d ' 



d 



d + c 1 —\f(c 1 )=I(f,f). 



(10) 



In the limit of small dissipation, the solution of fllOl) approaches a Maxwellian, i.e. /(( 



7T 



-d/2 



exp(— c 2 ). Therefore, a systematic approximation of the isotropic function f(c) 



can be found by expanding it in a set of Sonine polynomials, i.e. 

!oo 
i + E«A(c 2 
P =i 

which satisfy the orthogonality relations 



(11) 



dc0(c)5 p (c 2 )^(c 2 ) = 5 pp ^ 



(12) 



where 5 pp r is the Kronecker delta and M p a normalization constant. For general dimension- 
ality d, the first few Sonine polynomials are 



S Q (x) = 1 
Si(x) = —x + \d 
S 2 (x) = \x 2 -\{d + 2)x + \d{d + 2). 

The coefficients a„ are polynomial moments of the scaling function: 



(13) 
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dcS p (c 2 )f(c) = j^iSpic 2 )). 



(14) 



In particular a\ = (2/d)(Si(c 2 )) = 0, because the temperature definition (||) implies (c 2 ) 
\d. Moreover, a 2 is proportional to the fourth cumulant of the scaling form /(c), i.e. 



a-2 



d(d + 2) 



(c 4 )-Mrf + 2) =| (4)-3(clf 



(15) 



where we have used the relation, (c 4 ) = 3(c 4 )/[d(d + 2)], valid for any isotropic distribution 

he). 

To determine the coefficients a p we construct a set of equations for the moments 



((f) = / dcc*7(c), 



(16) 
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by multiplying (|T0| ) with (p = 1,2, . . .) and integrating over ci. For the moments fi p , 
defined in Eq. (|[), we obtain 




= f p(<?), (17) 

where the second line has been obtained by partial integration. For p = 2 the above equation 
reduces to a trivial identity because of the definition of temperature. 

The quantities /i 2 , \i p and (c p ) all depend on the unknown scaling function /(c). To calcu- 
late a 2 from (|T7j) we set p = 4, approximate the scaling form by /(c) = 0(c) {1 + a 2 5 , 2 (c 2 )}, 
and evaluate /z 2 , /Z4 and (c 4 ). The procedure is explained in more detail in the appendix 
and yields for general dimensionality d: 

16(1 - 2c 2 ) 

° 2 9 + 2Ad + 8ad- Ala + 30(1 -a)a 2 ' { ' 

This result for a 2 is plotted in Fig. p] as a function of a. 

In principle one can continue this approximation scheme by setting /(c) = 
4>(c) {1 + a 2 S' 2 (c 2 ) + ^^(c 2 )}, and then using (|i~7|) for p = 4 and p = 6 to obtain two 
coupled equations for a 2 and 03, and solve the resulting equations to obtain better approxi- 
mations for a 2 and a 3 than the previous ones, i.e. a 2 in (|l^) and = 0. As a 2 is already 
quite small, we do not calculate any higher coefficients a p (p > 3) in flTT|). 

For the three-dimensional case Goldshtein and Shapiro have calculated the coefficient 
a 2 and find the result 

qGS = 16(l-a)(l-2a 2 ) 
2 401 - 337a + 190(1 - a)a 2 ' 1 ' 

This result is only correct to linear order in 1 — a as the authors made an error in their 
algebraic calculations^. Their coefficient |a 2 | ~ 0.04 for all a G (0, 1), whereas the correct 
coefficient obeys |a 2 | ^ 0.2 for all a. However, the conclusion of Ref. || that the homoge- 
neous scaling form is well approximated by a Maxwellian remains valid for a large range of 
coefficients of restitution (say 0.6 ^ a < 1). For these values we have |a 2 | ^ 0.04 in three 
dimensions and |a 2 | ^ 0.024 in two dimensions. Our result for a 2 has been quantitatively 
confirmed by the Direct Simulation Monte Carlo results of Brey et al. 0. This will be 
discussed in section 0. 

To obtain the time dependence of the temperature, it is convenient to introduce the new 
time variable r representing the average number of collisions suffered per particle in a time 
t, and defined as dr = u (T(t))dt. This yields 

r = -ln(l + 7*/* ). (20) 
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In the unpublished appendices to their article, Eq. (E.10) should read A2 = 96 + 90a 2 . 
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Here to = l/^o(To) is the mean free time at the initial temperature T(0) = To. Next we find 
from Eq. (§) 



r(t ) = T exp(-2 7T ) = ITT ^. (21) 

In Eq. (|A.6| ) of the appendix, we derive for the cooling rate 7 = / dVL d )jj, 2 : 

7 = 7o{l + ^a 2 }, (22) 

where 70 = (1 — a 2 )/2d. Sela and Goldhirsch || have performed a numerical perturbation 
expansion of the Boltzmann equation to first order in e = 1 — a 2 and found the result 
7 = 7o (l - 0.0258e + C(e 2 )), which is close to the result 7 = 7„(1 - 3e/128 + C(e 2 )) = 
7o(l — 0.0234e + 0(e 2 )), obtained here. The method of appendix A also enables us to 
calculate the average collision frequency u = u[f] in the homogeneous scaling state with the 
result 



UJ = UJq 



{l (23) 



where the Enskog frequency ojq is defined in (^). Since the contribution from a 2 to 7 and u 
are small for all a, (0) and fl23"D are very well approximated by 70 and uq, respectively. 



III. UNIFORMLY HEATED SYSTEM 

To study this system we start from the stochastic equations of motion 

at m 

where Fj is the force due to collisions and $ >i is the random acceleration due to external 
forcing, which is assumed to be Gaussian white noise and uncorrelated for different particles, 
i.e. 

{L{t)i jP {t')) = ^SijScfi&it - 0, (25) 

where £g is the strength of the correlation. The validity of the above equations is based on 
the following assumptions: (i) the system is thermodynamically large, so that the condition 
J2i £i(t) = 0, imposed in computer simulations to guarantee momentum conservation in finite 
systems, can be ignored; (ii) the time between random kicks is small compared to the mean 
free time to, and therefore much smaller than the characteristic cooling time to/7 [ see Eq. 

The Enskog-Boltzmann equation for the single particle distribution function f(r,v,t) of 
a system heated in this way is corrected with a Fokker-Planck diffusion term (see e.g. Ref. 
T3fl), representing the change of the distribution function caused by the small random kicks, 



and reads in the spatially homogeneous case: 

2 Idvi 



dj(v u t) = X I(f, f) + f (■£-) /(vx, t). (26) 
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The diffusion coefficient £q is proportional to the rate of energy input |£q per unit mass. 
The equation for the temperature balance can be derived from Eq. ( [ZED m a similar fashion 
as in Eq. (||) for the cooling granular fluid, and reads 

^ = mil ~ 27^oT. (27) 



We are looking for a stationary solution of fl2"E|), where the heating exactly balances the 



loss of energy due to collisions, and the temperature becomes time independent. Again it is 
convenient to introduce a scaled distribution function by 

fM = -j(-), (28) 
where now the thermal velocity vq is time independent. Stationarity of / then requires 

T ^ + ^d^^kp^=°- (29) 

By multiplying this equation by cf and integrating over ci, we obtain the following set of 
equations which couple the moments (c p ~ 2 ) of the distribution to the moments \x v of the 
collision term, defined in Eq. (H): 



2 



r p(p + d-2)(<r 2 )=p p . (30) 

For p = 2 we recover the energy balance of Eq. (E7|), yielding for the stationary value of the 
thermal velocity in terms of fi2'- 



1/3 



Vo= —^TT • (31) 



Note that in order to obtain a finite temperature in the limit a — > 1, the a limit should be 
taken together with the limit £q ~~ ¥ 0- The above expression is used to write Eq. (^) in the 
form 

gp(p + d-2)(0=/V (32) 

Since a x = by definition of the temperature, i.e. (c 2 ) = |d, the first correction to Gaussian 
behavior is coming from ci2. To calculate it, we take p = 4 in Eq. (|3~2"P, use expression ( |A.8|) 
for /X4, and solve for 02 to finally obtain the result 

a 16(1 -a)(l -2a 2 ) 

2 73 + 56d -24ad- 105a + 30(1 - a)a 2 ' 1 J 

This function is shown in Fig. |2| for the two- and three-dimensional case. Again we find only 
small corrections to a Maxwellian distribution (02 < 0.086 in two dimensions and 0.067 in 
three). Therefore to a good approximation, \ii is given by its zeroth order approximation 
and the stationary temperature is found from Eqs. (|3T|) and ( |A.6|) as 

T ° = m (l-a^xW-. j ' (34) 
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IV. HIGH ENERGY TAILS 



In this section we will derive the asymptotic solution of the Enskog-Boltzmann equation 
(f29|) for high velocities in case the granular fluid is uniformly heated. Esipov and Poschel 
H have given a similar derivation for a freely evolving gas and found a high energy tail 
f(c) ~ exp(— Ac). The derivation in both cases proceeds along similar lines. If particle 1 
is a fast particle (ci ^> 1), the dominant contributions to the collision integral are collisions 
where particle 2 is typically in the thermal range, so that c 12 in the collision integral /(/, /) 
in (|9]) can be replaced by Ci. The gain term I g of the collision integral / can then be neglected 
with respect to the loss term Ii, as will be verified a posteriori at the end of this section. The 
collision integral /(/, /) then reduces to Ii —f3\Cif(ci), with f3\ = n^ 1 ^ 2 /T(\(d + 1)) as 
given in Eq. (|A.3|) of the appendix, and Eq. fllOD simplifies to 



j j(d + c^Kc) = -frcf{c). (35) 

The first term on the left hand side can be neglected with respect to the right hand side, 
and the large c solution has the form 

f(c)~AeM-—c), (36) 

where A is an undetermined integration constant. This solution corresponds to a tail which 
is overpopulated when compared to exp(— c 2 ). 

To determine the high energy tail of f(c) for the uniformly heated system, we proceed 
in a similar fashion and use (]3~ID to write Eq. (53) as 



™ » + 'M = - (37) 

For large velocities ci, the collision integral can again be replaced by — fliCif(ci), and Eq. 
( j37l) reduces to 

where we have used isotropy of the distribution function. Inserting solutions of the form 
f(c) oc exp(— Ac 3 ), we obtain the large c solution with B — | and A = §y^~, which is 
the only solution that vanishes for c —>■ oo. Again we find an enhanced population for high 
energies. 

To show that for c\ ^> 1 the gain term can be neglected with respect to the loss term, 
we use the asymptotic collision dynamics 

cT = c 1 -!(l + a- 1 )(c 1 -<x)<x 

c* 2 * = c 2 + |(l + a- 1 )(c 1 -6-)6-, (39) 

where we have replaced c 12 by Ci. If |cj • «r| 1, as is typically the case, c 2 in (^) can be 
neglected and we have 
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cT = c lv /l - 1(1 + a- x )(3 - a" 1 )^ ■ <x) 2 

^* = |(l + a- 1 )c 1 |ci-<r|»l, (40) 

where Cx is a unit vector. To demonstrate that f(c) ~ exp(— Ac B ) is a consistent large c 
solution, both in the freely evolving case with B = 1 and in the heated case with B = |, we 
compare the factor f(ci*)f(c%*) in J g with the factor f{c\)f [02) in 7; for large c, i.e. 

f(cT)f(c? 



f(ci)f(c 2 ) 



exp{-Ai(cr) B + (c?) B -cf}}. (41) 



The exponent is proportional to cf and strictly negative for a < 1 and B < 2, except 
for grazing collisions, where it vanishes. This happens inside a small 9 interval J of length 
0(l/ci) near = 7r/2, where |cr<x| = ci cos ~ Outside this interval the factor in fl4"I|) 

vanishes exponentially fast. Inside the interval J the factor in fl4T|) is 0(1). The contribution 
of this interval to the gain term can be estimated as Jj d6ci cos 0f{c\) ~ /(ci)/ci, where 
c x cos6> ~ 0(1)- Consequently, ~ 1/c 2 for large c\. 

In summary we have shown that /(c) ~ exp(— Ac B ) is a consistent large c solution of 
the Boltzmann Eqs. ( |10D and ( 29f) with B = 1 for the freely evolving fluid and 5 = | for 



the heated fluid. 



V. COMPARISON WITH SIMULATIONS 



In Refs. 1 11. 12], the undriven fluid of inelastic hard disks has been studied by molecular 



dynamics simulations. As long as the system is spatially homogeneous, measurements of 



the temperature decay confirm the validity of the homogeneous cooling law (21) where the 
cooling rate 7 is given by its zeroth order approximation 70 = e/2d. Also in the initial 
homogeneous state, the measured number of collisions C among N particles in a time t is 
consistent with r = 2C/N where r is given by Eq. (p0[), implying that the collision frequency 
io is very well approximated by its Enskog value u Q . 

So far, molecular dynamics simulations have not been able to obtain sufficient statistical 
accuracy to determine the fourth moment or the high energy tail of the velocity distribution. 
Such measurements are possible, however, by means of the Direct Simulation Monte Carlo 
method for the Enskog-Boltzmann equation. Using this method, Brey et al. |J have solved 
the nonlinear Boltzmann equation (|1|) for homogeneously cooling inelastic hard spheres 
(d = 3) and measured the fourth and sixth moment of the distribution f(c). Again the 
measured temperature decay shows no deviations of the cooling rate 7 from its Gaussian 
value 7q. Fig. 5 of Ref. compares their simulation data for the fourth cumulant 02 with 



first derived in |T3], and shows quantitative agreement. In particular, we predict that 



the fourth cumulant vanishes for a = I/a/2, which is very close to the value observed in 
the simulations. Also note that the simulation results disagree with the prediction of Ref. 
0. Moreover, the approximation /(c) = 0(c){l + c^S^c 2 )} shows a good agreement with 
the simulation data for the functional form of / (see Figs. 7 and 8 of Ref. [0]). This second 
Sonine approximation is qualitatively similar to the form presented in Fig. 3 of Ref. ||, 
calculated numerically to order 0(e). 
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It is also interesting to compare our theoretical predictions with recent molecular dynam- 
ics results of Peng and Ohta || on the heated granular fluid. These authors have measured 
the temperature relaxation T(t) in a fluid of N inelastic hard disks of mass m = 2 at an 
area fraction = ^na 2 N/L 2 ~ 0.16 and heating rate £ 2 = (5V) 2 /3t h ~ 1.67 x 10~ 4 , 
where the randomly added velocity components are sampled from a uniform distribution 
on the interval (-SV, 5V) where 5V = 10~ 3 , measured in system length L per unit time, 
and th = 2 x 10~ 3 is the period between random kicks. For the pair distribution of hard 



disks at contact, X-> we use the approximate form |~5| x = (1 — -^0)/(l — (p) 2 ~ 1.32. The 



16 ■ 

-3 



steady state temperature predicted by Eq. (|34J) then becomes T = 1.15 x 10~ 3 for a = 0.8 
The simulations yield Tq 1111 ~ 1.21 x 10~ 3 (see Fig. 1 of Ref. ||), in fair agreement with the 
Enskog theory. 

Moreover, Eq. ( |34"D predicts that T depends on the heating rate £g = (5V) 2 /3th and 
inelasticity e = 1 — a 2 , as 

The measurements show an exponent A = 0.65 ± 0.01. The theoretical prediction ([34]) gives 
Co — 0.092. The simulation result of Ref. M Cq° ~ 5.0 x 10 -3 , corrected^ with the conversion 
factor (L/er) 2 / 3 , gives <f im ~ 5.0 x 10~ 3 x (70.9) 2//3 ~ 0.086, again in fair agreement with the 
theoretical prediction. 

Eq. (|27|) also gives a prediction for the approach of T{t) to T . By observing that tu oc \/T 
and linearizing Eq. (|27| ) around T , one obtains the solution T{t) = T +Ti exp(— t/r ) with 
tq = 2To/3m£ 2 . For the parameters <fi = 0.16 and a = 0.8 of Fig. 1 in Ref. this yields 
ro = 2.3, and their simulations yield TQ im = 1/b' = 2.6. 

Next, we compare the collision frequency ujq in (0) from the Enskog theory for elastic hard 
disks with the collision frequency u; sim , measured in Ref. f8|. If there are C binary collisions 
among iV particles in a time t, then the collision frequency is 2C/Nt. The simulation 
results at a = 0.2, 0.4, 0.6, 0.7 are respectively u sim ~ 2.93, 1.67, 1.49, 1.49, and the 
Enskog predictions for the same a values are ujq ~ 1.17, 1.22, 1.33, 1.44. The Enskog 
frequency ojq ~ V^o decreases according to ([42]) with increasing e = 1 — a 2 , whereas u; sim 
increases strongly with e. The simulation results for u sim suggest that the Enskog theory 
gives reasonable predictions for a > 0.6. Similar conclusions have been obtained by Orza et 



al. |L2[ for the homogeneous cooling state of a freely evolving fluid of inelastic hard disks. 
The deviations at larger inelasticities are probably caused by clustering and the onset of 
kinetic collapse, which strongly increases the collision frequency. 

Finally, we observe that the overpopulation in the high energy tail ~ exp(— Ac 3 / 2 ) of the 
steady state distribution function has also been observed in the simulations of Ref. 0, but 
their statistical accuracy is too low to make any quantitative comparison. 



2 The parametrization ( |42D of the measurements in Ref. || with Cq° ~ 5.0 x 10~ 3 does not 
reproduce their data points at a = 0.8. Moreover, the corresponding £ value ~ 1.4 x 10~ 3 ) is 
far away from any data point included in their Fig. 6. 
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VI. CONCLUSIONS AND OUTLOOK 



We have investigated non-Gaussian behavior in granular fluids of smooth inelastic hard 
spheres or disks, both in the absence of an external forcing and in a system uniformly heated 
by random accelerations. In a freely evolving granular fluid, we find for all inelasticities very 
small corrections to the cooling rate and the collision frequency due to non-Gaussian char- 
acteristics of the homogeneous cooling state. As a consequence, such deviations have never 
been observed in computer simulations on homogeneous systems. Our result for the fourth 
cumulant in the homogeneous cooling state has been confirmed by the computer simula- 
tions of Brey et al. |]T|. These authors used the Direct Simulation Monte Carlo method 
to obtain an accurate homogeneous solution of the nonlinear Enskog-Boltzmann equation. 
This method might make it feasible to obtain quantitative information on the high energy 
tail ~ exp(— Ac) to test the theoretical predictions. It certainly could be used to measure 
the fourth moment of the homogeneous solution in a uniformly heated system, a quantity 
calculated in the present paper. Again, in this case we predict very small corrections to 
the stationary temperature and collision frequency due to non-Gaussian properties of the 
homogeneous stationary state, which are possibly too small to measure in computer simula- 
tions. Moreover, it would be interesting to investigate the validity of our prediction for the 
overpopulation of the high energy tail ~ exp(— Ac 3 / 2 ) with A = \\J^j^- 

Long range spatial correlations, measured in one- and two-dimensional simulations of 
heated granular fluids J7|-[J, are currently being analyzed using the ring kinetic equations 
corresponding to kinetic equation (|26"D with the Fokker-Planck diffusion term, as well as 
by fluctuating hydrodynamics with external noise. The approach of adding external noise 
has some similarity with the Edwards- Wilkinson model Jl6|| for the growth of a granular 



surface on which particles are impinging at random. As a consequence, long range spa- 
tial correlations in the velocity-velocity and density-density correlation functions are to be 
expected. 

The clustering reported in the simulations of Ref. ||, which causes an enhancement of 
the collision frequency with respect to the Enskog value for higher inelasticities (a < 0.6), 
has not yet been explained in terms of a linear or nonlinear stability analysis of the long 
wavelength hydrodynamic modes of the system. 

T.v.N. acknowledges support of the foundation 'Fundamenteel Onderzoek der Materie 
(FOM)', which is financially supported by the Dutch National Science Foundation (NWO). 

VII. APPENDIX A 

In this appendix we calculate the quantities /i 2 ,/i 4 and (c 4 ), which are required in (|17D 
to calculate the coefficient a 2 in ( |TTD by setting /(c) = 0(c) {1 + g^S^c 2 )}, where 0(c) is 



the Maxwellian. In fact, the moment (c 4 ) in ([IB]) requires only moments of the Gaussian 
distribution. A straightforward calculation gives 

(c 4 ) = J dcc 4 0(c){l + a 2 S 2 (c 2 )} 

= \d(d + 2){l + a 2 }. (A.l) 
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Next we consider the moments /i p (p = 2,4) in (§). With the help of Eq. (|3]) it can be 
transformed into 



(A.2) 



= ~\ I dci / dc 2 f d<r(ci2 • <x)0(ci)0(c 2 ) x 
{l + a 2 [S 2 (cl) + S 2 (c 2 2 )} + 0{al)} A(c? + eg), 

where the operator A is defined below Eq. ([3]). In the following, terms of 0(a 2 ) will be 
neglected. For a > 0.6 this is a safe approximation as can be checked from the results Eqs. 
and PD for a 2 . 



To evaluate (|A.2| ) we introduce center of mass and relative velocites by Ci = C + |c 12 
and c 2 = C — hc 12 . Moreover, we need the angular integral 



d<r(ci 2 ■ or) 

■n/2 



„ J"d<X(cos0) n 

P 

J d<T 



in £_ de ( sin6, ) d ^ 2 ( cos9 ) n tir(2±i) 



(A.3) 



2 °" fl j^Mtfitae)*-' " ' r(s±*) : 

where c 12 = c 12 /|c 12 | is a unit vector. Using the relations between Gaussian moments, it is 
straightforward to derive the relations: 

(c 12 ■ *C 2 AC n c? 2 ) = \{d + n)(c 12 ■ <xAC™c™> 

(ci 2 • <xC 4 AC n c™> = ±(d + n)(d + n + 2)(c 12 -&AC n c? 2 ) , (A.4) 

where 

(V(c 12) C))o = J dc 12 ^-^exp(-ic? 2 )/dC(^) d/2 exp(-2C 2 )^(c 12 ,C) (A.5) 

denotes a Gaussian average over c 12 and C. The above formulas are very helpful in evaluating 
the moments \i p in (|A.2 ). With the help of ( |A.3| ) and ( |A.4j ) one finds 

to = K 1 -a 2 )^(c? 2 )o {1 + ^2} 

= 5(i-« 2 )^r{i + ^}- ( A - 6 ) 

To calculate // 4 we need the quantity 

A(c 4 + c 4 ) = 2(1 + a) 2 (c 12 ■ <x) 2 (C ■ <x) 2 + ±(a 2 - l) 2 (c 12 • <r) 4 



+( a 2 -l)(c 12 -<x) 2 C7 2 + ±i 



a 



1)(C12 • 0-) 2 C 2 2 



-4(l + a)(c 12 -o-)(C-o-)(C-c 12 ) 
One finds after long and tedious calculations 

A*4 = /^3(c? 2 )o {7i + a 2 T 2 } , 

with 

T 1 = i(l-a 2 )(d+| + a 2 ) 



(A.7) 
(A.8) 

(A.9) 



T 2 = ^(1 - a 2 )(10d+ 39 + 10a 2 ) + f (1 + a)(d - 1). 

For the homogeneous cooling solution, inserting the results ( |A.1|) , ( |A.6j ) and (|A.8|) into (|l 
for p = 4 yields a closed equation for a 2 . Neglecting again small contributions 0(a 2 ), we 
solve for a 2 , and the result in flT8|) is recovered. Eq. (|33|) corresponding to uniform heating 
is found by inserting ( |A.6|) and (|A.8|) into ( |32"1) for p = 4. 
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FIG. 1. Fourth cumulant a-i versus a for homogeneous cooling solution in a freely evolving fluid. 
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FIG. 2. Fourth cumulant a,2 versus a for the stationary state of a uniformly heated system. 
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